Maria Doyle, Jessica Chung, Vicky Perreau, Kim-Anh LĂȘ
Cao
09 October 2023
R for Biologists course
R takes time to learn, like a spoken language. No one can expect to
be an R expert after learning R for a few hours. This course has been
designed to introduce biologists to R, showing some basics, and also
some powerful things R can do (things that would be more difficult to do
with Excel). The aim is to give beginners the confidence to continue
learning R, so the focus here is on tidyverse and visualisation of
biological data, as we believe this is a productive and engaging way to
start learning R. After this short introduction you could use this
book to dive a bit deeper.
Most R programmers do not remember all the command lines we share in
this document. R is a language that is continuously evolving. They use
Google extensively to use many new tricks. Do not hesitate to do the
same!
There are two ways to work in RStudio in the console or in a script.
We can type a command in the console and press Enter to run
it.
Exercise: Run the command below in the console.
1 + 1
[1] 2
Or we can use an R script. To run a command line, we place the cursor
on the line you want to run, and then either:
Click on the run button on top of the panel
Use Ctrl + Enter (Windows/Linux) or
Cmd + Enter (MacOS).
We can also highlight multiple lines at once and run them at
once.
Exercise: Create a script, from the top menu in
RStudio: File > New File > R Script, then type the
command below in the script and run it with one or the two options.
2 + 2
[1] 4
â Using a script is recommended practice
If we work in the console, we do not have a good record (recipe) of what
we have done. It is good practice to use an R script to be able to
re-use the code and remember what we have done in the past.
Commenting
Comments are notes to ourself or others about the commands in the
script. they are useful also when you share code with oters. Comments
start with a # which tells R not to run them as
commands.
# testing R
2 + 2
[1] 4
Keeping an accurate record of how youâve manipulated your data is
important for reproducible
research. Writing detailed comments and documenting your work are
useful reminders to your future self (and anyone else reading your
scripts) on what your code does.
In the main exercise we will use Rmarkdown to create
reproducible and beautiful analysis reports.
Working directory
Opening an RStudio session launches it from a specific location. This
is the âworking directoryâ. R looks in the working directory by
default to read in data and save files. You can find out what the
working directory is by using the command getwd(). This
shows you the path to your working directory in the console. It will be
in this format: /path/to/working/directory (Mac) or
C:\path\to\working\directory (Windows). We recommend your
data and R scripts in the same directory.
To define the working directory, in RStudio click in the menu on
Session > Set Working Directory > Choose directory
and choose your folder.
Exercise: Create a folder for this course somewhere
on your computer. Name the folder for example,
Intro_R_course. Then, to set this folder as your working
directory. Set your working directory as indicated above, then click
âOpenâ. The files tab will now show the contents of your new folder.
Exercise: Save the script you created in the
previous section as intro.R in this directory. You can do
this by clicking on File > Save and the default location
should be the current working directory
(e.g. Intro_R_course).
R Packages or libraries
A package is a collection of functions along with documentation,
code, tests and example data created by R developers who wish to share
their methods and code to others. Often these packages are not installed
by default on your computer.
We use the following command line to install a package:
install.packages("tidyverse")
We then load the package in our working directory:
library("tidyverse")
Packages in the CRAN or Bioconductor
Packages hosted on the CRAN
(stands for Comprehensive R Archive Network) are often generic package
for all sorts of data and analysis. Bioconductor is an ecosystem that
hosts packages specifically dedicated to biological data. The
installation of packages is a bit different, e.g to install the
mixOmics package we type:
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("mixOmics")
You dont need to remember this command line, as it is featured in the
Bioconductor package page (see here
for example).
One advantage of Bioconductor packages is that they are well
documented, updated and maintained every six months.
Here is how the example vignette looks like:
browseVignettes("mixOmics")
R functions
In mathematics, a function defines a relation between inputs and
output. In R (and coding languages) it is the same. A function takes
several inputs called arguments inside parentheses, and output
some results.
Exercise: Compare these two outputs. In the second
line we use the function sum.
2+2
[1] 4
sum(2,2)
[1] 4
Operation assignment
It is useful to store data or result so that we can use them later on
for other parts of the analysis. We can use either the operator
= or <-. In both cases the object where we
store the result is on the left-hand-side, and the result from the
operation is on the right-hand-side.
Exercise: Compare the two outputs.
result1 = 2+2
result1
[1] 4
result2 <- 2+3
result2
[1] 5
â Nomenclature of objects
We recommend you use explicit naming conventions of your objects, for
example data.raw1 and data.normalised rather
than data1 and data2 so that you can remember
various steps in your analysis.
I need help!
To see what any function in R does, type a ? before the
name and help information will appear in the Help panel on the right in
RStudio.
?sum
What is really important is to scroll down the examples to understand
how the function can be used in practice. You can use this command line
to run the examples:
example(sum)
Google and Stack Overflow
are also useful resources for getting help.
Common R errors
R error messages are common and often cryptic. You most likely will
encounter at least one error message during this tutorial. Some common
reasons for errors are:
Case sensitivity. In R, as in other programming languages, case
sensitivity is important. ?install.packages is different to
?Install.packages.
Missing commas
Mismatched parentheses or brackets or unclosed parentheses,
bracketsand apostrophes
Not quoting file paths
When a command line is unfinished, the â+â in the console will
indicate it is awaiting further instructions. Press ESC to cancel the
command.
To see examples of some R error messages with explanations see here
Letâs get started with data!
In this tutorial, we will learn some R through creating plots to
visualise data from an RNA-seq experiment.
The GREIN
platform (GEO RNA-seq Experiments Interactive Navigator) provides
>6,500 published datasets from GEO that have been uniformly
processed. It is available at http://www.ilincs.org/apps/grein/. You can search for a
dataset of interest using the GEO code. GREIN provide QC metrics for the
RNA-seq datasets and both raw and normalized counts. We will use the
normalized counts here. These are the counts of reads for each gene for
each sample normalized for differences in sequencing depth and
composition bias. Generally, the higher the number of counts the more
the gene is expressed.
RNA-seq dataset form Fu et al.
Here we will create some plots using RNA-seq data from the paper by
Fu et
al. 2015, GEO code GSE60450. This study examined expression in basal
and luminal cells from mice at different stages (virgin, pregnant and
lactating). There are 2 samples per group and 6 groups, 12 samples in
total.
Tidyverse
The tidyverse package is a collection of R packages that
includes the extremely widely used ggplot2
package.
The tidyverse makes data
science faster, easier and more fun.
Load the package
We use library() to load in the packages that we need.
As described in the cooking analogy in the first screenshot,
install.packages() is like buying a saucepan,
library() is taking it out of the cupboard to use it.
Exercise: Load your tidyverse library. If you get an
error message, it means that you have not installed it! (see the code in
the previous section).
library(tidyverse)
Load the data
The files we will use are csv comma-separated, so we will use the
read_csv() function from the tidyverse. There is also a
read_tsv() function for tab-separated values.
Exercise: Download the data.zip file here. Unzip
the file and store this data folder in your working
directory.
We will use the counts file called
GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv stored
in the data folder. The path to the file should be
data/GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv.
Exercise: Load the data into your R working
directory. We will store the contents of the counts file in an
object called counts.
# read in counts file
counts <- read_csv("data/GSE60450_GeneLevel_Normalized(CPM.and.TMM)_data.csv")
New names:Rows: 23735 Columns: 14ââ Column specification âââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââ
Delimiter: ","
chr (2): ...1, gene_symbol
dbl (12): GSM1480291, GSM1480292, GSM1480293, GSM1480294, GSM1480295, GSM1480296, GSM1480297, GSM1480298, GSM1480299,...
âč Use `spec()` to retrieve the full column specification for this data.
âč Specify the column types or set `show_col_types = FALSE` to quiet this message.
# read in metadata
sampleinfo <- read_csv("data/GSE60450_filtered_metadata.csv")
New names:Rows: 12 Columns: 4ââ Column specification âââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââââ
Delimiter: ","
chr (4): ...1, characteristics, immunophenotype, developmental stage
âč Use `spec()` to retrieve the full column specification for this data.
âč Specify the column types or set `show_col_types = FALSE` to quiet this message.
Do not be afraid of the outputs. There is some information output by
read_csv on âcolumn specificationâ It tells us that there is a missing
column name in the header and it has been filled with the name ââŠ1â,
which is how read_csv handles missing column names by default. It also
tells us what data types read_csv is detecting in each column. Columns
with text characters have been detected (col_character) and also columns
with numbers (col_double). We wonât get into the details of R data types
in this tutorial but they are important to know when you get more
proficient in R. You can read more about them in the R
for Data Science book.
Looking at the data
When assigning a value to an object, R does not print the value. We
do not see what is in counts or sampleinfo.
But there are ways we can look at the data.
Exercise: Option 1: Click on the counts
object in your global environment panel on the right-hand-side of
RStudio. It will open a new Tab.
Exercise: Option 2: type the name of the object and
this will print the first few lines and some information, such as number
of rows.
sampleinfo
We can also use dim() to see the dimensions of an
object, the number of rows and columns.
We can also take a look the first few lines with head().
This shows us the first 6 lines.
head(sampleinfo)
We can look at the last few lines with tail(). This
shows us the last 6 lines. This can be useful to check the bottom of the
file, that it looks ok.
tail(sampleinfo)
Dimension of the data
It is very useful to check that what you have read in R corresponds
to the number of lines and columns of the data;
Exercise: Check how many rows and columns are in
sampleinfo.
dim(sampleinfo)
[1] 12 4
This show us there are 12 rows and 4 columns.
In the Environment Tab in the top right panel in RStudio we can also
see the number of rows and columns in the objects we have in our
session.
Column and row names of the data
Exercise: Check the column and row names used in in
sampleinfo. Comment on the results you get.
Subsetting is very important in R to extract parts of the data you
want to analyse. In R, a data set has the form of a matrix with rows and
columns. You can access these rows, or columns by adding brackets.
For example, here we extract the first row â1â of the data, as
indicated on the left-hand-side of the colon:
sampleinfo[1,]
Here we extract the firsecond column â2â of the data, as indicated on
the right-hand-side of the colon:
sampleinfo[,2]
You can use a combination of number of row and column to extract one
element in the matrix:
sampleinfo[1,2]
Exercise: Work out what the following commands mean,
then include those to subset sampleinfo (either on the rows
or the columns):
1:3
[1] 1 2 3
c(1, 3)
[1] 1 3
â The c() function
We use this function extremely often in R when we have multiple items
that we are combining (âcâ stands for concatening). We will see
it again in this tutorial.
In the previous section, when we loaded in the data from the csv
file, we noticed that the first column had a missing column name and by
default, read_csv function assigned a name of â...1â to it. We can
change this column to something more descriptive by combining what what
we have learnt.
Exercise: Describe what the next command line is
doing:
Exercise: Do the same to rename the first column
name from â...1â to âgene_idâ in counts.
colnames(counts)[1] <- "gene_id"
Other useful commands: str and
summary
str() shows us the structure of our data. It shows us
what columns there are, the first few entries, and what data type they
are e.g. character or numbers (double or integer).
summary() generates summary statistics of our data. For
numeric columns (columns of type double or integer) it outputs
statistics such as the min, max, mean and median. We will demonstrate
this with the counts file as it contains numeric data. For character
columns it shows us the length (how many rows).
summary(counts)
gene_id gene_symbol GSM1480291 GSM1480292 GSM1480293 GSM1480294
Length:23735 Length:23735 Min. : 0.000 Min. : 0.000 Min. : 0.00 Min. : 0.00
Class :character Class :character 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.00 1st Qu.: 0.00
Mode :character Mode :character Median : 1.745 Median : 1.891 Median : 0.92 Median : 0.89
Mean : 42.132 Mean : 42.132 Mean : 42.13 Mean : 42.13
3rd Qu.: 29.840 3rd Qu.: 29.604 3rd Qu.: 21.91 3rd Qu.: 19.92
Max. :12525.066 Max. :12416.211 Max. :49191.15 Max. :55692.09
GSM1480295 GSM1480296 GSM1480297 GSM1480298 GSM1480299 GSM1480300
Min. : 0.00 Min. : 0.00 Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Median : 0.58 Median : 0.54 Median : 2.158 Median : 2.254 Median : 1.854 Median : 1.816
Mean : 42.13 Mean : 42.13 Mean : 42.132 Mean : 42.132 Mean : 42.132 Mean : 42.132
3rd Qu.: 12.27 3rd Qu.: 12.28 3rd Qu.: 27.414 3rd Qu.: 26.450 3rd Qu.: 24.860 3rd Qu.: 23.443
Max. :111850.87 Max. :108726.08 Max. :10489.311 Max. :10662.486 Max. :15194.048 Max. :17434.935
GSM1480301 GSM1480302
Min. : 0.000 Min. : 0.000
1st Qu.: 0.000 1st Qu.: 0.000
Median : 1.629 Median : 1.749
Mean : 42.132 Mean : 42.132
3rd Qu.: 23.443 3rd Qu.: 24.818
Max. :19152.728 Max. :15997.193
â Multiple methods
There are multiple ways to do things in R. For example we used the
function colnames in the prevous exercise, but we could
have used the rename() function. When programming, you will
often find many ways to do the same thing. Often there is one obvious
method depending on the context.
Formatting data
Converting from wide to long format
We will first convert the data from wide format into long format to
make it easier to work with and plot with the ggplot
package.
Instead of multiple columns with counts for each sample, we just want
one column containing all the expression values, as shown below:
We can use pivot_longer() to easily change the format
into long format.
Explanation: We use cols = starts_with("GSM")
to tell the function we want to reformat the columns whose names start
with âGSMâ. pivot_longer() will then reformat the specified
columns into two new columns, which weâre naming âSampleâ and âCountâ.
The names_to = "Sample" specifies that we want the new
column containing the columns we sepcified with cols to be
named âSampleâ, and the values_to = "Count" specifies that
we want the new column contining the values to be named âCountâ.
We could also specify a column range to reformat. The command below
would give us the same result as the previous command.
Alternatively, we could specify the columns we do not want
to reformat and pivot_longer() will reformat all the other
columns. To do that we put a minus sign â-â in front of
the column names that we donât want to reformat. This is a pretty common
way to use pivot_longer() as sometimes it is easier to
exclude columns we donât want than include columns we do. The command
below would give us the same result as the previous command.
Exercise: Type each command line above, then look at
the data. Check their dimension (dim()) and length
(length())
seqdata
Joining two tables
Now that weâve got just one column containing sample ids in both our
counts and metadata objects we can join them together using the sample
ids. This will make it easier to identify the categories for each sample
(e.g. if itâs basal cell type) and to use that information in our
plots.
Explanation: We use the function full_join()
and give as arguments the two tables we want to join. We add
by = join_by(Sample == sample_id) to say we want to join on
the column called âSampleâ in the first table (seqdata) and
the column called âsample_idâ in the second table
(sampleinfo) when the values match:
allinfo <- full_join(seqdata, sampleinfo, by = join_by(Sample == sample_id))
Exercise: Examine the code above, run and then look
at the data.
allinfo
Exercise
View the help page of the head function and find the
âArgumentsâ heading. What does the n argument do? How many
rows and columns do you get with
head(sampleinfo, n = 8)?
Store the output of the first 20 lines of the counts
object in a new variable named subset_counts. What is the
gene_symbol in row 20?
View the values in the GSM1480291 column from your
subset_counts object using the $ subsetting
method. What is the 5th value?
View the help page of the mean function. What is the
mean of the column of values you obtained from the previous
question?
How can you use pivot_longer to transform
dat into a âtidyâ data called dat_long that
contains 3 columns (sample, experiment, count).
Join dat_long to sample_info using the
common column called âsampleâ. How many rows do you get if you use i)
full_join, ii) left_join, iii)
right_join, iv) inner_join?
Answers
The n argument specifies the number of elements the
head() function will return.
head(sampleinfo, n = 8) will return the first 8 rows of the
sampleinfo object.
subset_counts <- head(counts, n = 20)
subset_counts$GSM1480291
137.8874 (which is obtained by running
mean(subset_counts$GSM1480291))
110 from running
nrow(full_join(dat_long, sample_info, by = join_by(sample == sample)))
20 from running
nrow(left_join(dat_long, sample_info, by = join_by(sample == sample)))
110 from running
nrow(right_join(dat_long, sample_info, by = join_by(sample == sample)))
20 from running
nrow(inner_join(dat_long, sample_info, by = join_by(sample == sample)))
Plotting with ggplot2
ggplot2 is a plotting package that
makes it simple to create complex plots. One really great advantage
compared to classic R packages is that we only need to make minimal
changes if the underlying data change or if we decide to change our plot
type, for example, from a box plot to a violin plot. This helps in
creating publication quality plots with minimal amounts of adjustments
and tweaking.
ggplot2 likes data in the âlongâ
format, i.e., a column for every variable, and a row for every
observation, similar to what we created with pivot_longer()
above. Well-structured data will save you lots of time when making
figures with ggplot2.
**â Ggplot architecture*
Ggplot graphics are built step by step by adding new elements using the
+. Adding layers in this fashion allows for extensive
flexibility and customization of plots.
To build a ggplot, we use the following basic template that can be
used for different types of plots. Three things are required for a
ggplot:
The data
The columns in the data we want to map to visual properties (called
aesthetics or aes) e.g. the columns for x values, y values
and colours
The type of plot (the geom_)
There are different geoms we can use to create different types of
plot e.g. geom_line()geom_point(),
geom_boxplot(). To see the geoms available take a look at
the ggplot2 help or the handy ggplot2
cheatsheet. Or if you type âgeomâ in RStudio, RStudio will show you
the different types of geoms you can use.
Creating a boxplot
Letâs plot boxplots to visualise the distribution of the counts for
each sample. This helps us to compare the samples and check if any look
unusual.
Note: with ggplot the â+â must go at the end of the line, it
canât go at the beginning.
Exercise: Run the following command line. Identify
the key functions aes and type of plot:
This plot looks a bit weird. Itâs because we have some genes with
extremely high counts. To make it easier to visualise the distributions
we usually plot the logarithm of RNA-seq counts. Weâll plot the Sample
on the X axis and log~ 2~ Counts on the y axis. We can log the Counts
within the aes(). The sample labels are also overlapping
each other, we will show how to fix this later.
We get a warning here about rows containing non-finite values being
removed. This is because some of the genes have a count of zero in the
samples and a log of zero is undefined. We can add +1 to every count to
avoid the zeros being dropped (âpsuedo-countâ).
The box plots show that the distributions of the samples are not
identical but they are not very different.
Exercise: violin plot
Box plots are useful summaries, but hide the shape of the
distribution. For example, if the distribution is bimodal, we would not
see it in a boxplot. An alternative to the boxplot is the violin
plot, where the shape (of the density of points) is drawn.
Letâs choose a different geom to do another type of plot. Using the
same data (same x and y values), try editing the code above to make a
violin plot using the geom_violin() function.
Answers
# Plotting a violin plot
ggplot(data = allinfo, mapping = aes(x = Sample, y = log2(Count + 1))) +
geom_violin()
Colouring by categories
Letâs add a different colour for each sample.
Exercise: Get the help file for
geom_boxplot and scroll down to Aesthetics
heading. It specifies that there is an option for colour.
In the code below we specify we want to map the Sample column to
colour =. As we are mapping colour to a column in our data
we need to put this inside the aes().
Colouring the edges wasnât quite what we had in mind. Look at the
help for geom_boxplot to see what other aesthetic we could
use. Letâs try fill = instead.
ggplot(data = allinfo, mapping = aes(x = Sample, y = log2(Count + 1), fill = Sample)) +
geom_boxplot()
That looks better. fill = is used to
fill in areas in ggplot2 plots, whereas
colour = is used to colour lines and points.
A really nice feature about ggplot is that we can easily colour by
another variable by simply changing the column we give to
fill =.
Exercise
Modify the plot above. Colour by other variables (columns) in the
metadata file:
characteristics
immunophenotype
`developmental stage` (As there is a space in the column name we
need to use backticks around the name (``).
Note: backticks are not single quotes (ââ). The backtick key is
usually at the top left corner of a laptop keyboard under the ESC key.
Check what happens if you donât use backticks.)
The geom_boxplot function can also take in additional
arguments. For example, you can decrease the size of the outlier points
by using the outlier.size argument like so:
geom_boxplot(outlier.size = 0.5). View the help page for
geom_boxplot. Can you find a way to hide outliers
altogether? Plot a boxplot with hidden outliers.
Answers
ggplot(data = allinfo,
mapping = aes(x = Sample, y = log2(Count + 1), fill = characteristics)) +
geom_boxplot(outlier.shape = NA)
Creating subplots for each gene using faceting
With ggplot we can easily make subplots using faceting. For
example we can make stripcharts.
These are a type of scatterplot and are useful when there are a small
number of samples (when there are not too many points to visualise).
Here we will make stripcharts plotting expression by the groups (basal
virgin, basal pregnant, basal lactating, luminal virgin, luminal
pregnant, luminal lactating) for each gene.
Shorter category names
First we use the function mutate() to add a column with
shorter group names to use in the plot, as the group names in the
characteristics column are quite long.
Exercise: Run the code below, then look at the data
using the head() function. What haschanged?
Exercise: Run the code above and look at
mygenes_counts, what did it do?
To identify these 8 genes, we used pipes (%>%)
to string a series of function calls together (which is beyond the scope
of this tutorial, but totally worth learning about independently!).
However, the boxplots are not very informative as we only have two
values per group. Letâs plot the individual points instead using the
geom_point() to make a scatterplot.
The points are overlapping so we will make a jitter plot using
geom_jitter() which adds a small amount of random variation
to the location of each point so they do not overlap. If is also quite
common to combine jitter plots with other types of plot, for example, jitter
with boxplot.
We can save plots interactively by clicking Export in the Plots
window and saving as e.g. âmyplot.pdfâ. Or we can output plots to pdf
using pdf() followed by dev.off(). We put our
plot code after the call to pdf() and before closing the
plot device with dev.off().
R markdown is a great tool embedded in RStudio that will produce
reproducible and very neat reports.
R Markdown is a specific type of file format (.Rmd)
designed to produce documents that include both code and text. We use it
in place of the R script, as you will see in the exercises
below.
â Folder organisation with R Markdown
The working directory should look like this:
a folder where we store the data
a folder where the figures from the report will be saved
the report in .Rmd (this is what we will work with)
the rendered report in .pdf
Exercise: Together, letâs go through the different
steps of compiling the .Rmd file to produce a .pdf report.
A good place to start with RMarkdown is with this cheatsheet.
Main exercises
Exercise 1: plotting
Here we will work on the RNA-seq data from Fu et al. However, feel
free to upload your own data instead and apply your new learning!
The raw counts for the RNA-seq data are available in the data folder
with the name âGSE60450_GeneLevel_Raw_data.csvâ.
Plot a boxplot. Do the samples look any different to the normalised
counts?
Make subplots for the same set of 8 genes. Do they look any
different to the normalised counts?
The normalised counts for another data set
âGSE63310_GeneLvel_Normalized(CPM.an.TMM)_data.csvâ are also available.
Plot boxplots and violin plots colouring the samples using different
columns in the metadata file.
# read in counts file
counts_GSE63310 <- read_csv("data/GSE63310_GeneLevel_Normalized(CPM.and.TMM)_data.csv")
colnames(counts_GSE63310)[1] <- "gene_id"
# read in metadata
sampleinfo_GSE63310 <- read_csv("data/GSE63310_filtered_metadata.csv")
colnames(sampleinfo_GSE63310)[1] <- "sample_id"
seqdata_GSE63310 <- pivot_longer(counts_GSE63310, cols = -c("gene_id", "gene_symbol"),
names_to = "Sample", values_to = "Count")
allinfo_GSE63310 <- full_join(seqdata_GSE63310, sampleinfo_GSE63310, by = join_by(Sample == sample_id))
# Colouring by categories
# characteristics
ggplot(data = allinfo_GSE63310, mapping = aes(x = Sample, y = log2(Count + 1), fill = characteristics)) +
geom_boxplot()
# strain/background
ggplot(data = allinfo_GSE63310, mapping = aes(x = Sample, y = log2(Count + 1), fill = `strain/background`)) +
geom_boxplot()
# age
ggplot(data = allinfo_GSE63310, mapping = aes(x = Sample, y = log2(Count + 1), fill = age)) +
geom_boxplot()
Exercise 2: reading data, cleaning and reformatting
This exercise is useful to apply your learnings on how to wrangle
data (a step that is often overlooked but takes a long time!)
We will look at the âbreast_tumors.csvâ, âgene_description.txtâ
and âgene_name.xlsxâ data files, then read each file using one of the
appropriate function read_csv(), read_tsv() or
read_excel() (from readxl). Name each file appropriately.
NOTE: you will still need to load readxl explicitly, because it is not a
core tidyverse package loaded via library(tidyverse).
Extract the treatment variable from breast_tumor.csv, inspect and
solve issues in this file (hint: there are some labelling
inconsistencies).
Combine the data from âgene_name.xlsxâ and
âgene_description.txtâ
Answers
# read the files
gene.exp<-read_csv("data/breast_tumors.csv")
gene.description<-read_tsv("data/gene_description.txt",skip = 2, col_names = FALSE)
colnames(gene.description)<-c("gene.ID", "description")
gene.name<-read_excel("data/gene_name.xlsx")
# extract the treatment variable
treatment<-as.factor(toupper(gene.exp$treatment))
# combine the gene.name and gene.description
gene.info<-full_join(gene.name, gene.description, by = join_by(gene.ID))
Session Info
The last thing weâll do run the sessionInfo() function.
This function prints out details about your working environment such as
the version of R youâre running, loaded packages, and package versions.
Printing out sessionInfo() at the end of your analysis is
good practice as it helps with reproducibility in the future.
sessionInfo()
Key Points
Tabular data can be loaded into R with the tidyverse functions
read_csv() and read_tsv()
Tidyverse functions such as pivot_longer(),
mutate(), filter(), select(),
full_join() can be used to manipulate data
A ggplot has 3 components: data (dataset), mapping (columns to plot)
and geom (type of plot). Different types of plots include
geom_point(), geom_jitter(),
geom_line(), geom_boxplot(),
geom_violin().
facet_wrap() can be used to make subplots of the
data
The aesthetics of a ggplot can be modified, such as colouring by
different columns in the dataset, adding labels or changing the
background
We might want to change the colours. To see what colour names are
available you can type colours(). There is also an R
colours cheatsheet that shows what the colours look like.
There are built-in colour palettes that can be handy to use, where
the sets of colours are predefined. scale_colour_brewer()
is a popular one (there is also scale_fill_brewer()). You
can take a look at the help for scale_colour_brewer() to
see what palettes are available. The R
colours cheatsheet also shows what the colours of the palettes look
like. Thereâs one called âDark2â, letâs have a look at that.
We can change the axis labels and add a title with
labs(). To change the x axis label we use
labs(x = "New name"). To change the y axis label we use
labs(y = "New name") or we can change them all at the same
time.
ggplot(data = mygenes_counts,
mapping = aes(x = Group, y = log2(Count + 1), colour = Group)) +
geom_jitter() +
facet_wrap(~ gene_symbol) +
labs(x = "Cell type and stage", y = "Count", title = "Mammary gland RNA-seq data")
Themes
We can adjust the text on the x axis (the group labels) by turning
them 90 degrees so we can read the labels better. To do this we modify
the ggplot theme. Themes are the non-data parts of the plot.
ggplot(data = mygenes_counts,
mapping = aes(x = Group, y = log2(Count + 1), colour = Group)) +
geom_jitter() +
facet_wrap(~ gene_symbol) +
labs(x = "Cell type and stage", y = "Count", title = "Mammary gland RNA-seq data") +
theme(axis.text.x = element_text(angle = 90))
We can remove the grey background and grid lines.
There are also a lot of built-in themes. Letâs have a look at a
couple of the more widely used themes. The default ggplot theme is
theme_grey().
ggplot(data = mygenes_counts,
mapping = aes(x = Group, y = log2(Count + 1), colour = Group)) +
geom_jitter() +
facet_wrap(~ gene_symbol) +
labs(x = "Cell type and stage", y = "Count", title = "Mammary gland RNA-seq data") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90))
ggplot(data = mygenes_counts,
mapping = aes(x = Group, y = log2(Count + 1), colour = Group)) +
geom_jitter() +
facet_wrap(~ gene_symbol) +
labs(x = "Cell type and stage", y = "Count", title = "Mammary gland RNA-seq data") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90))
There are many themes available, you can see some in the R graph
gallery.
We can also modify parts of the theme individually. We can remove the
grey background and grid lines with the code below.
The groups have been plotted in alphabetical order on the x axis and
in the legend (that is the default order), however, we may want to
change the order. We may prefer to plot the groups in order of stage,
for example, basal virgin, basal pregnant, basal lactate, luminal
virgin, luminal pregnant, luminal lactate.
First letâs make an object with the group order that we want.
Next we need to make a column with the groups into an R data type
called a factor. Factors in R are a special data type
used to specify categories, you can read more about them in the R for Data
Science book. The names of the categories are called the factor
levels.
Weâll add another column called âGroup_fâ where weâll make the Group
column into a factor and specify what order we want the levels of the
factor.
Take a look at the data. As the table is quite wide we can use
select() to select just the columns we want to view.
select(mygenes_counts, gene_id, Group, Group_f)
Notice that the Group column has <chr> under the
heading, that indicates is a character data type, while the Group_f
column has <fct> under the heading, indicating it is
a factor data type. The str() command that we saw
previously is useful to check the data types in objects.
str(mygenes_counts)
str() shows us Group_f column is a Factor with 6 levels
(categories).
We can check the factor levels of a column as below.
levels(mygenes_counts$Group_f)
The levels are in the order that we want, so we can now change our
plot to use the âGroup_fâ column instead of Group column (change
x = and colour =).
We could do similar if we wanted to have the genes in the facets in a
different order. For example, we could add another column called
âgene_symbol_fâ where we make the gene_symbol column into a factor,
specifying the order of the levels.
Exercise
Make a colourblind-friendly plot using the colourblind-friendly
palettes here.
Create a plot (any plot whatsoever) and share it with the class
by pasting the image in the Google Docs link provided in your workshop.
You plot should use the subtitle argument in the
labs function to add a unique identifier (e.g. a message
and your name or initials) which is displayed below the title.
Tip: An easy way to copy your plot in RStudio is using the plot
paneâs export option and selecting âCopy to ClipboardâŠâ. You can then
paste it into the provided Google document.
Answers
# Define colourblind colours
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2")
# An example of a colourblind-friendly stripchart
ggplot(data = mygenes_counts,
mapping = aes(x = Group_f, y = log2(Count + 1), colour = Group_f)) +
geom_jitter() +
facet_wrap(~ gene_symbol) +
labs(x = "Cell type and stage", y = "Count", title = "Mammary gland RNA-seq data") +
theme(axis.text.x = element_text(angle = 90)) +
theme_bw() +
scale_colour_manual(values = cbPalette)
# An example of a colourblind-friendly boxplot
ggplot(data = allinfo,
mapping = aes(x = Sample, y = log2(Count + 1), fill = characteristics)) +
geom_boxplot() +
theme_bw() +
scale_fill_manual(values = cbPalette)
